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FOREWORD 


This is a final report on the research project, "Numerical Solutions of 
Three-Dimensional Navier-Stokes Equations for Closed-Bluff Bodies," for the 
period January 1, 1987 to December 31, 1987. Special attention during this 
period was directed to "Topology and Grid Adaption for High-Speed Flow 
Computations." The work was supported by the NASA Langley Research Center 
(Computer Applications Branch of the Analysis and Computation Division) 
through the cooperation agreement NCC1-68. The cooperative agreement was 
monitored by Dr. Robert E. Smith, Jr. of ACD-Computer Applications Branch. 
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TOPOLOGY AMD GRID ADAPTION FOR HIGH-SPEED FLOW COMPUTATIONS 


By 


J. 


S. Abo! hassani 1 and S. N. Tiwari^ 


SUMMARY 

This study investigates the effects of grid topology and grid adaption on 
numerical solutions of the Navier-Stokes equations. In the first part of this 
study, a general procedure is presented for computation of high-speed flow 
over complex three-dimensional configurations. This includes the grid 
generation and solution algorithm for Navier-Stokes equations in a general 
three-dimensional curvilinear coordinate system. The flow field is simulated 
on the surface of a Butler wing in a uniform stream. Results are presented 
for Mach number 3.5 and a Reynolds number of 2,000,000. The 0-type and H-type 
qrids have been used for this study, and the results are compared together and 
with other theoretical and experimental results.^ The results demonstrate that 
while the H-type grid is suitable for the leading and trailing edges, a more 
accurate solution can be obtained for the middle part of the wing with an 0- 
type grid. In spite of some discrepancies, the present numerical results 
compare favorably with the experimental results. In the second part of this 
study, methods of grid adaption are reviewed and a method is developed with 
the capability of adapting to several variables. This method is based on a 
variational approach and is an algebraic method. Also, the method has been 
formulated in such a way that there is no need for any matrix inversion. This 
method is used in conjunction with the calculation of hypersonic flow over a 
blunt-nose body. A movie has been produced which shows simultaneously the 
transient behavior of the solution and the grid adaption. 

For both cases, the simulations are done by integrating the viscous 
Navier-Stokes equations. These equations govern the unsteady, viscous, 
compressible and heat-conducting flow of an ideal gas, and all viscous terms 
are retained. The equations are written in curvilinear coordinates so that 
the body surface is represented accurately. The computer codes are written in 
FORTRAN, is vectorized and currently run on the CDC Vector Processing System 
(VPS-32, CYBER 205) computer. The results indicate the viability and validity 
of the proposed methods. 


^Research Associate, Dept, of Mechanical Engineering and Mechanics 
Eminent Professor, Dept, of Mechanical Engineering and Mechanics 
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Chapter 1 


INTRODUCTION 

Continuum problems in engineering are quite often modeled by 
systems of nonlinear partial differential equations. These equations 
are usually complex and, in most instances, must be solved by numerical 
means. The numerical solution of the governing partial differential 

equations has two steps: (1) grid generation and (2) numerical 

integration. Grid generation is the division of the solution domain 
into discrete interconnected points called grid points. The accumula- 
tion of grid points is called a grid. The location of grid points is 
primarily a function of boundary geometry and the physics of the 

problem. The grid should conform to the boundaries and be concentrated 
in regions where there are large gradients. In the second step, the 
derivatives in the partial differential equations are approximated with 
algebraic expressions (usually by Taylor series expansions where higher 
order terms are truncated). The most commonly used techniques for 

integrating the governing equations are classified as finite difference 
techniques, finite volume techniques or finite element techniques. The 
accuracy of a numerical solution depends on both the solution technique 
and the grid. Grid points must be appropriately defined to apply 
boundary conditions and must be sufficiently close together to resolve 
the physics to the desired level of accuracy. They must also be 

oriented relative to each other in such a manner that errors are not 



introduced into the solution. On the other hand, computer speed and 
memory limit the number of grid points that can be used in the solution 
of a given problem. It is, therefore, necessary to distribute grid 
points to maximize overall accuracy while covering the entire region of 
interest. 

If the grid points are ordered in such a way that the relationship 
between any grid point and its neighbors is the same for all grid 
points, then this grid is called a structured grid. A structural grid 
can be generated numerically by determining the values of the physical 
coordinates in the physical domain from the values on the boundaries. 
This can be accomplished in two basic ways: (1) by algebraic 

interpolation from the boundary values, and (2) by solving a set of 
partial differential equations with the boundary geometry as a boundary 
condition. The resulting boundary-fitted coordinate system is a 
curvilinear coordinate system having some coordinate lines (surfaces in 
three-dimensions) conforming to the shape of each boundary. When the 
governing equations are transformed onto such a coordinate system, a 
finite approximation can be made using neighboring points at coordinate 
line intersections, without the need for interpolation, regardless of 
boundary shape and boundary movement. Thus, quite general code can be 
written for the numerical solutions of the governing partial 
differential equations on arbitrary regions. Although many of the 
accomplishments in grid generation have occurred within the field of 
computational fluid dynamics, the techniques are equally applicable in 
electromagnetics, solid mechanics, and other areas involving solutions 
of partial differential equations in an arbitrary region. The 
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literature on grid generation is extensive and this is critically 
reviewed in [1, 2]*. 

An examination of the Taylor series expansion of a function about a 
point in the solution domain, reveals that the truncation error depends 
on derivatives of the solution and characteristics of the grid, such as 
distribution of spacing, orthogonality and aspect ratio. The selection 
of a grid topology has direct effect on the solution of a given problem 
through the introduction of singularities and constraints on the 
orthogonality of the grid. Eriksson [3] has studied the effects of grid 
singularities on the solution of the Euler equations. It was pointed 
out that the nonconservative centered scheme is likely to be unstable at 
mesh singularities, whereas, the conservative centered finite volume 
scheme is stable in a local sense. Also, it was concluded that the 

; i 

local time-stepping gives rise to exponentially growing modes for 
nonconservative schemes. Mastin and Thompson [1, 4] examined two 

sources of truncation error in the numerical solutions of partial 
differential equations on a curvilinear coordinate system. The error 
sources are derived from grid spacing and the degree of nonorthogonality 
(grid skewness). It is possible for a poor distribution or orientation 
of grid points to Introduce errors into a numerical solution. For 
example, sudden changes in the line spacing and excessively skewed lines 
can Introduce negative numerical diffusion into a solution. Although 
precise orthogonality is not essential, some error terms vanish for an 
orthogonal system. Also, Raithby [5] has observed that an excessive 


* 
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grid skewness will exaggerate the truncation error. The literature 
survey indicates a lack of information on the effects of grid topology 
on the solution, especially in three-dimensions. In general, the 
coordinate system should have lines concentrated in regions of an 
expected high variation of the physical solutions. The coordinate 
system should be coupled with the physical solution so that the 
coordinate lines continually adapt to resolve the evolving gradients in 
the physical solution. A significant amount of work on the grid 
adaption is available in the literature and this is discussed, in 
detail, in Chap. 6. 

In regard to the grid generation and solution of continuum 
problems, this study has two distinct objectives. These are: (1) the 

qualitative assessment of errors resulting from the relative orientation 
of grid points and overall grid topology, and (2) the proper 
redistribution of the grid points over a region to minimize the 
truncation error through grid adaption. In order to study the effects 
of grid topology, flow over a Butler wing [6] has been simulated using 
two different grid topologies. Solutions obtained with the two 
topologies are compared with each other and with experimental results 
obtained by Squire [7-9]. In the second part of this study, a method of 

grid adaption with the capability of adapting the grid points to several 

variables is proposed. This method is formulated in such a way that it 
is not necessary to solve a system of equations and is, therefore, very 
efficient computationally. 

The physical models used in this study are discussed in Chap. 2. 

The formulations of the governing equations are presented in Chap. 3. 

The method of solution and the initial and boundary conditions are 
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explained briefly in Chaps. 4 and 5, respectively. Chapter 6 contains a 
review of the grid-adaption method and discussion of the proposed new 
method. Finally, some critical results are presented in Chap. 7. 
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Chapter 2 


PHYSICAL MODELS 

The two basic physical models considered for this study are 
discussed briefly in this chapter. These are a Butler wing and a blunt 
leading edge of a panel holder. 

The Butler wing is a good test case for investigating the effects 
of grid topology on the numerical solutions of Navier-Stokes equations. 
This is due to a unique feature of its geometry. The Butler wing is a 

f 

delta wing which was proposed by D. S. Butler [6]. The plan form of the 
wing is an isosceles triangle, and the leading edges of the wing lay 
along the Mach lines of the unperturbed stream. The first twenty 
percent of the wing is conical and the last eighty percent of the wing 
has elliptical cross sections with increasing eccentricity along the x- 
axis (Fig. 2.1). At the trailing edge, the elliptical cross section has 
infinite eccentricity and is a straight line. The Butler wing is 
symmetric about (x-z) and (x-y) planes. This permits the use of one 
quarter of the entire physical domain with zero angle of attack (Fig. 
2.1). However, if the angle of attack is greater than zero, then half 
of the physical domain should be considered. The semi-major and minor 
axes are given by 
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1.1 
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Fig. <^.la Physical model of a Butler Wing (top view) 





Major axis 

(semi -span) = j 


0 < x < L 

(2.1a) 

Minor axis 

(thickness on = j 


0 < x < 0.2L 

(2.1b) 


centerline) i 





-io- 

rX - 0 • 2L 1 4 ] 
L ~0.8x J J 

0.2L < x < L 

(2.1c) 


where 



The model is 0.8 ft. (0.2438 m) long, and the geometry has been 
generated for a Mach number of 3.5. That is the semi-apex angle of the 
plan form and the initial conical nose is sin -1 (l./3.5) = 16.602°. 

Butler has compared the experimental results for surface pressure 
with the theoretical results [6]. These theoretical results are 
obtained from inviscid equations of motion which are simplified by using 
the slender-body approximations. Walkden and Caine estimated the 
pressure on the surface of a Butler wing at zero incident in a steady 
uniform stream. They numerically integrated the two semi -character is tic 
forms of the equations governing Inviscid supersonic flow of an ideal 
gas with constant specific heats [10]. Squire has obtained experimental 
results for a Butler wing with varying Mach number and angle of attack 
[9], In all previous analytical and numerical investigations, the 
Inviscid form of the equations of motion has been used. 

Grid generation is the first step which should be considered in 
obtaining flow field solutions over any configuration. Due to the data 
base management of the present program, it is necessary to map an entire 
physical domain into a rectangular parallelepiped. Among the grid 
types, selection of an 0-type grid for cross sections in the stream-wise 
direction produces a point singularity at the nose tip and a line 
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singularity at the trailing edge (Fig. 2.2). Nevertheless, an O-type 
grid maps the entire solid boundary onto an entire face of the 
parallelepiped. It is also possible to generate orthogonal grid in the 
regions where there is relative high curvature. However, an H-type grid 
does not map the solid boundary onto an entire face of the computational 
box (Fig. 2.2). This creates a potential problem in updating the 
boundary conditions near the leading edges of the wing and also the grid 
in some regions could be highly skewed. But, there are no singularities 
in the grid. Figure 2.2 shows a topological comparison between H-type 
and 0-type grids for the Butler wing. Eoth types of grids have been 
used in this study, and the results are compared with other numerical, 
analytical ana experimental results. 

Another aspect of this study is the grid adaption for high-speed 
flow computation. In the hypersonic flow about blunt bodies (Fig. 2.3), 
the temperature, pressure and density of the flow increase almost 
explosively across a shock wave. At the same time, the curved shock 
wave is close to the body. Numerical simulation of this phenomena has 
been a great challenge to the computational fluid dynamics researchers. 
Presently, there is a great deal of interest in improving the quality of 
numerical simulation techniques, and grid adaption is one way to achieve 
thi s goal . 

The accuracy of finite-difference solutions depends on the fineness 
of the grid. Therefore, the finer the grid, the more accurate the 

numerical solution will be. Also, the accuracy of solutions depends on 
the resolution of the solution gradient. The presence of large 
gradients causes the error to be large in the difference approximation 
of derivatives. In the presence of shock waves, more artificial 
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edge of a panel holder 


diffusion must be added to retain adequate smoothness of the solutions. 
Therefore, there is a need for schemes that can resolve large gradients 
without adding additional grid points. An adaptive scheme moves the 
grid points to regions of high gradients, when the locations of these 
gradients are not known a priori. Also, an adaptive method reduces the 
total number of grid points required to achieve a given accuracy, but it 
requires more computer time. In some instances, the computer time makes 
grid adaption impractical. The ideas used in the construction of 
adaptive grid techniques are limited only by one's imagination; and any 
scheme that works in the sense of providing a better solution is a good 
one. The ultimate answer to numerical solutions of partial differential 
equations may well be to dynamically adapt grids, rather than to devise 
more elaborate difference representations and solution methods [11]. 

In both cases, the flow is simulated by solving the Navier-Stokes 
equations numerically. The equations are unsteady, compressible, 

viscous and three- or two-dimensional. The time dependency of the 
governing equations allows the solution to progress naturally from an 
arbitrary initial guess to an asymptotic steady state, if one exists. 
The equations are transformed from physical coordinates to computational 
coordinates, allowing the solutions to be computed in a rectangular 
domain. The equations are solved by the MacCormack time-split technique 
[12, 13] which is vectorized and programmed to run on the CDC VPS-32 
(CYBER 205) computer. The code is written in 32-bit (half-word) 

FORTRAN. The details of the formulation and solution procedure are 
presented in the subsequent chapters. 
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Chapter 3 


GOVERNING EQUATIONS 


3.1 Navier-Stokes Equations 

The governing equations for a thermal fluid system are the conser- 
vation of mass, momentum and energy. These equations are developed for 
an arbitrary region under the assumption that the system is a continuum. 
Equations of motion for a viscous, compressible, unsteady and heat- 
conducting fluid can be written as [14] 


Continui ty : 

3 p 

Tt + 

V • 

(pu) =0 , 

(3.1a) 

Momentum: 

&(pu) A 

bt 

V • 

( puu - x) = 0 , 

(3.1b) 

Energy: 

3(E) , 

IT + 

V • 

(Eu + q - u • x) = 0 , 

(3.1c) 


where E is the total energy per unit volume given by E = p(e + vv/2) 
and e is the internal energy per unit volume. Equations ( 3 . lb-3 . lc ) can 
be simplified by assuming that the stress at a point is linearly 
dependent on the rate of strain (deformation) of the fluid (Newtonian 
fluid) 


5u . &u • 3u. 

T ij = ' P 6 ij + ^ + + 6 ij * 


(3.2a) 


The Kronecker delta function is denoted by 6.., and p' is the second 

* J 

coefficient of viscosity which is related to the coefficient of bulk 
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viscosity (k) by the expression k = 2p/3+ti'. The contribution of < can 
be neglected If the pressure in a fluid is not changed abruptly during 
its expansion or contraction, in other words, the hydrostatic pressure 
is assumed to be equal to the average of the normal stresses. Under 
this assumption, the stress tensor can be related to the pressure and 
velocity components as 


'u * • p 6 u + “ 


5u. 

J 


du. 


2 

1 


1J 


du ki 


(3.2b) 


For an isotropic system, the heat flux in Eq. (3.1c) can be expressed in 
terms of temperature gradient (Fourier law of heat conduction) as 

q = - K 7 T (3.3) 


where K is the coefficient of thermal conduction. A common 
approximation used for viscosity is based on the kinetic theory of gases 
usino an idealized intermoleculai force potential. The relation is 


where 




3/2 T o * S o 


S . 198.6 °R 
0 

T = 492 °R 
0 

u = 0.35xl0 6 (lbf-sec)/ft 2 


(3.4) 


the coefficients of thermal conduction K can be determined from the 
Prandtl number as 
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(3.5) 


yuC 


where C y Is the specific heat at a constant volume and y is the ratio 
of specific heats. 

It is essential to have a supplementary relation to close the 
system of equations, Eqs. (3.1a-3.1c). By neglecting the intermolecular 
forces (thermally perfect system), the thermodynamic properties can be 
described by the equation of state 

P = pRT (3.6) 

where R is the universal gas constant. The assumption of thermally 
perfect gas permits the internal energy to be expressed as a function of 
temperature only i.e., e=e(T). In addition, the assumption of a 

) * f J 

calorically perfect gas [e(0)=0] allows the following relation 

e = C y T . (3.7) 

A substitution of Eq. (3.6) into Eq. (3.7) results in 


P = pe (y-1) . 


(3.8) 


The equations of motion are in conservative form. For simplicity, these 
equations can be written in a compact vector form as 


where 


au 

3T 





(3.9) 
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X 

1 

pvv 
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yy 

pvw 

- T 

yz 

Ev + 

V 


+ p 


4 ^ + PV, 
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pw 


puw 

x zx 

pvw 

T zy 

pWW 

' T zz 

£w + 

• 

q z - 


+ P 


xx 


= -P + 2 n - | 

K 3x 3 ^ v ^x oy TEr » 


_ ,du . dv. 

xy " ^ %y + * 


_ _ ,,/3w , dti. 

V - “<37 + 37> • 


yy 


-P + 2a ^ .,(® u + j. ® w \ 

P ^ 3 y 


yz 


,dv , dw> 


x - _p + 9 m _ 2 / 5 u , dv 5 Wi 

zz P 2 “ 5 z 1 ** ( - 5 x + 3 y + T 5 z> - 


a * u \x + "xy + WT xz ' 


$ = UT + VT + Wx 

y xy yy w yz * 


*z = UT xz + VT yz + WT zz * 


« = - K 
q x K dx * 


o = -K dT 
q y K W * 


o - -K 
q z K "5z • 
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3.2 Navier- Stokes Equations in Computational Coordinates 


For the sake of generality, the governing equations are transformed 
from a physical domain into a computational domain as 



+ 



9F 3G 5H\ 

"5C ~5Zf 


= 0 . 


(3.10) 


The stress tensor, dissipation function and heat conduction must also be 
transformed from a physical domain into a computational domain. Using 
Eqs. (A. 23c) and (A. 24a), they can be expressed as 


x 


xx 


- _p + 2 U. &U 3n 3u 3C 9u-| /-3C 3u 3n 3u 

" p z ^ 3 x - 5 Z 37 33 + 37 2/3 l 37 3 £ + -57 33 


ac au 
37 ac 


, ac av an av ac av . ac aw an aw 

+ 373c 3y 33 37 3c 35-£c + -5z-53 


ac aw-, 

3z 3T j 


9 


T 


xy 


/-ac au an au ac au ac av . an av 

' ** l 37 33 + 37 33 37 3£ + 37 3f + 37 33 


+ 


ac av^ 
37 W 


* 


T 


xz 


/•ac aw an aw ac aw ac au 
* l 37 3£ 3733 373£ + 3z3^ 


an au ac au-i 

3z 33 3z 3£J 


9 


yy 


-P + 2p ( 


ac av 
37 3£ 


, an av 
37 33 


. ac av 
37 3? 


) - 2/3 


r ac au 
l 37 3 ? 


an au 

37 33 


ac au 
"57 ‘Sc 


ac av an av ac av ac aw an aw 

3y "5? "57 33 "57 "5c 3z ac az an 


ac aw^ 
az ac J 


» 
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yz 


= n ( 


3C av 
3z Sc 


an dv 

37 ^n 


ac &V . dC dw 

<5“z ‘W, "5y "5f 


an dw ac dw\ 
+ TSy oy 
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= -P + 2n ( 


ac aw 

"5z ^ 


a-n aw 


. ac aw 
+ ISz ^ 


) - 2/3 |1 ( 


ac au 

"5x "5c 


, at) au 

~5x^i 


ac au 
ax ac 


ac av 5 t) av ac av 

+ Ty "5y 3n 3y ^ 


ac aw an aw 

+ -51 *5? '5n 


ac aw 

“5z ST 


) 


9 


i, f ac aT 
= " K ^TSc 


3ti 3T 

"Sx "5n 


+ 


ac aT> 

•5x TZ ] 


- -K ( 


ac aT . ar) aT 
3y K W ^ 


ac ai\ 

TyW 


9 


- -K ( 


ac aT 

Tz M 


an aT 
+ -55 ^ 


+ 


ac dT-i 

SzW * 


The transformation is based on the chain rule. The transformation 
coefficients can be computed from a functional relation between the 
computational coordinates and the physical coordinates as 

C = C (x, y, z) , 

n = n (x, y, z) , (3.11) 

C = C (x, y, z) . 

If the relation in Eq. (3.11) were known, the transformation coefficient 

could have been computed by direct differentiation. If not, after some 
algebraic manipulations (Eq. A. 14), the transformation coefficients can 
be computed by 
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ac 

ac 

ac ) 

( ax 

dx 

ax 

37 

3y 

Bz / 


■Bn 

BC 

an 

an 

an \ 

= CJ] = 

ay 

ay 

37 

■By 

Bz / 

an 

"BC 

ac 

ac 

ac l 

1 az 

az 

az 

37 

3y 

Bz ) 

f al 

"Bn 

BC 


-1 


(3.12) 


where [J - *] is defined as 


[J] « 


rdy bz 
l 3B TC " 

(by bz 

" l 3B Be ~ 

r by bz 
[ bt dn " 


ay bz 
3 S ’Bn 

3y bz 

K K 


) -(■ 


5y bZ} 

Bn 3c j 


3x bz 

Bn Be 


3x bz 

B£ Bn 


) 


z' dx bz dx bz 

l 3c 3C ” Be B£ 

r ax bz ax 
^ac an ~ an Be 


\ r ax ay ax ay^, 

J l B3 "5c " B£ W 

^ r ax ay ax ay-v 

] “ l 3B Be " Be 3c J 

i _ &x ayi 

' (• dc an ^ * 


Bn ac 

(3.13) 


IJ _1 I ■ 


ax ax ax 
ac 'Bn ac 

' ay ay ay 
\3f "Bn "Bt/ 

' bz bz bz\ 
13£ "Bn BIT , 


ax r 3y bz by az-v 

~ 31 1 an Be ~ Bt Bn J 

, ax /-ay az ay az\ 
+ "5? l Bc 'Bn "Bn B £ ’ 


ax ,-ay az ay az^ 

"Bn '■Bt Be Bt ac 


In the case of a Butler wing, the grid planes are perpendicular to 
the x-coordinate. Consequently, physical coordinates can be written as 

x = x (C) , 

y = y (C, n, C) , (3.14) 

z = z (c n, C) . 

This reduces the transformation coefficients from nine to five non-zero 
elements and, therefore, reduces the memory requirements. 
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In the case of grid adaption, the physical geometry is two- 
dimensional. Therefore, all derivatives with respect to the x and l 
coordinates are set to zero. This simplifies the governing equations to 
the following form 


where 


5U 

Ft 


+ [\ ( 

T1 2 J 


dG 5H 
Fn Fn 


) + [ C c y ] ( 

g z 


5G dH-i _ n 

FC Fc J ~ 



G = 


pv 

PW - T yy + P 

pVW - T yz 

E V + q y - 4> y + Pv 


(3.15) 


H 


pw 

PVW - T zy 

pww - 1 Z1 
Ew + q z - 


+ P 

4 > + Pw 

z 


There are four non-zero transformation coefficients. 

In both cases, the transformed governing equations are called Chain 
Ruled Conservation Law Form (CRCLF) [15]. However, the governing 
equations can be written with metric coefficients inside the 

differentiations; this form is called the Strong Conservation Law Form 
(SCLW) . It has been shown [15] that CRCLF requires no special 

considerations on how to compute the metric coefficients or their 
derivatives and it also requires that fewer arithmetic operations be 
performed compared to the other forms. Also, it has been shown that it 
has the ability to capture weak shocks. 
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Chapter 4 


METHOD OF SOLUTION 


A time-marching method is used to compute the solution so that the 
possible transient features can be readily captured. This explicit 
method is a time-split predictor-corrector algorithm which is second- 
order accurate in time and space [12]. The governing equations are 
split into three groups of operators, each aligned with transformed 
coordinates. Then, these equations, Eq. (3.10), are discretized in the 
computational directions. In a compact form, they can be expressed as 


,n+l 


i » j.k = [ L t/ A VJ f L c (At C^ [ L ^ A V] [ L c ( A t c )] [L^ (At^ ) ] 


where 


At n = At c = 7 


(4.1) 


and L^, , and L c are the operators in the 5 , n, and c directions, 

respecti vely . A time step is completed in this algorithm with the 
application of each operator applied symmetrically about the middle 
operator. Operator can be defined as 


for the predictor step: 

At, 


L (At ) = U 0Ut 
V V i,j,k ’ 


(4.2a) 


U 


i.j.k = U i"j,k ‘ TT [(F i ‘ F i- 1 ) H i + (G i " G i- 1 } § 1 


+ 


(H. - H 



as 


i] 


j.k 


(4.2b) 
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for the corrector step: 




+ {G i+l " G i } *5y 1 + (H i+l " V H 1 \ J 

J » K 


Operator can be defined as 

L (At) = U? u1 : . , 
r) n i > J * k 


for the predictor step: 

At 

0. . = U 1 ." . . - -r-2- [ (F . 

i.j.k i.j.k St) L j 


F . ,) ^2 j + (G . 
j-1 J ' j 


r i 5ti i 

G j-r dy J 


J J 
for the corrector step: 


+ (H .i " W S rt uk • 


..out = 1 r yin y . J 1 f(F 

U i,j,k ?' U i,j,k i,j,k At) j +1 


At 


F ) On i 
j J " 5 x J 


+ (G j+l " V W j + ( Vl ‘ V *51 


Operator can be defined as 


W ■ u ?j. k • 


for the predictor step: 


. At- 

u 4 . . = u: n . . - T J ’ [ ( f. 
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for the corrector step: 
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+ (G k + l - G k>W k + (H k+l - H k>i k ]) • (4 - 4c » 

' <3 

Fluxes are computed by a forward-difference approximation for the 
predictor and a backward-difference approximation for the corrector. 
Therefore, the algorithm is second-order accurate In space. More 
details can be found in [12-14]. 


The solution is stable if the time step of each operator does not 
exceed the allowable step size for that operator. The finite-difference 
scheme is consistent if the sums of the time steps for each operator are 
equal. The solution is second-order accurate if the operators are 
applied symmetrically [12-14], 


This method has a time-step stability limit, but there is no 
rigorous stability analysis available. A commonly used conservative 
time- step is 


« < mtn [M ♦ lli + |iJ. + c (_ij._lj._lj)] 


-1 


(4.5) 


Ax Ay~ A z' 

where c is the local speed of sound. This is valid for cartesian 
coordinates. Using Fourier series, a similar equation is derived for 
general curvilinear coordinates in Appendix B 
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In the supersonic and hypersonic regions, there exists large 
gradients which require a very fine grid to resolve them. Most central- 
difference methods admit a solution which has sawtooth or plus-minus 
waves with the shortest wave length that the grid can support. In the 
case of a nonlinear problem, these short waves interact, vanish, and 
reappear again as distorted long waves, or oscillations. These 
oscillations eventually cause the solution to blow up if they are not 
resolved. These numerical oscillations are caused by truncation error 
and can be reduced by grid refinements. The oscillations of "low 
frequency" can be suppressed by adding a fourth-order damping term. A 
common damping used is pressure damping. This is usually expressed in 
cartesian coordinates as 




The last term in Eq. (4.7a) is computed by a forward-difference 

approximation for the predictor step and a backward-difference 
approximation for the corrector step. It is common to use only an 
approximation of the second derivatives of pressure in the transformed 
coordinate system, instead of using the actual rigorous transforma- 
tions. The use of covariant velocities may raise problems for complex 

grid topologies. One remedy for general curvilinear coordinates is to 
use contravariant velocities defined by 

222 1/2 

U = (£ x u + £ y v + s z w) (x^ + + z ? ) 

V = (n x u + n y v + ti z w) (X* + yjj + (4.7b) 
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The equations for g^, g^ and g^ are defined in Eq. (A. 10), and U, V 
and W are defined in Eq. (4.7b) 
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Chapter 5 


INITIAL AND BOUNDARY CONDITIONS 

In computational fluid dynamics the Initial conditions usually 
correspond to a real situation for a transient problem, or a rough guess 
for a steady state problem. In practice, initial conditions are 
obtained from experiments, empirical relations, approximate theories or 
previous computational results. An Improper initial guess may result in 
generating unrealistic strong transient waves which propagate throughout 
the computational region, dominating the flow field and eventually 
leading to a solution failure. An important requirement for the initial 
conditions is that they should be physically as close as possible to the 
actual nature of the flow field in the region under study. This will 
minimize the number of iterations required for convergence. An 
attractive approach Is to initialize the entire flow field with a crude 
and simple guess (e.g., free stream conoition). During the course of 
the computation, both body and upstream boundary conditions are changed 
in a gradual manner to their final values over a prescribed number of 
iterations. In the present study, this technique Is applied in only one 
step which is equivalent to Impulsive initial conditions. 

It is equally important to implement a realistic, accurate and 
stable method to determine boundary conditions. The application of 
certain conditions may cause numerical instability even though the flow 
is physically stable. Most of the boundary conditions currently 
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Implemented are drawn mainly from Intuition, simple analytical 
expressions, wind tunnel experiments and computational experimentation. 
In the selection of boundary conditions, consideration should be given 
to the following criteria: convergence, stability, computer time and 
above all the physical justification. For the Butler wing case, there 
are five different boundary conditions. They are upstream, downstream, 
outer, solid boundaries, and symmetry. 

5.1 Upstream Boundary Condition for a Butler Wing 

For the case of an H-grid, the upstream boundary is located at six 
grid point spacings ahead of the nose of the wing. The following 
undisturbed free stream conditions are assumed for this boundary 

Vo.k ■ • (5 - 1 > 

For the case of an 0-type grid, the upstream boundary is set at 
five percent of the chord from the tip of the wing to avoid the point 
singularity. The conical assumption has been made for this boundary 
[16]. Flow is said to be conical if the physical conditions such as 

pressure and velocity do not vary with position along any ray through a 
point, referred to as the vertex. For this case, the viscous-conical 

solutions are obtained for a cone at the proper angle of attack. This 
is done by creating a conical grid which has straight lines (rays) from 
the vertex (Fig. 5.1). Then, the conical Navier-Stokes equations are 
Integrated for the middle plane (plane B), then planes A and C are set 
equal to plane B. This procedure is repeated until convergence is 

reached. This solution is valid provided the wing is sufficiently 

slender. 
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c. 



Flg - 5.1 Conical grid 
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5.2 Downstream Boundary Conditions for a Butler Wing 

A zero gradient in the ^-direction (parallel to the primary 
direction of flow) is assumed for the downstream boundary, i.e.. 


au 


= o . 

IL.j.k 


(5.2a) 


A backward-difference is used to approximate the Eq. 
results in 

V,j,k = V-l,j,k * 


(4.2a) which 
(5.2b) 


5.3 Far-field Boundary Conditions for a Butler Wing 

The outer boundary is located far away from the body to avoid any 
influence on the interaction region. Presently, a zero normal gradient 
of fluxes is assumed for this boundary, i.e.. 


5U 

Sn 


= 0 

i ,JL,k 


(5.3a) 


Similar to the down stream boundary condition, a backward difference is 
used to approximate Eq. (5.3a) which results in 


U i,JL,k " U i,JL-l,k 


(5.3b) 


5.4 Solid-Wall Boundary Conditions for a Butler Wing 

The walls are assumed to be impermeable and no-slip boundary 
conditions are applied, therefore, all velocity components are assumed 
to be zero. Similarly, the wall is assumed to have a constant 
temperature T w . A zero normal pressure gradient is assumed for the 
solid surface, i.e. 
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&p 

3n 


= 0 . 


(5.4) 


1.1. k 


This appears to be a boundary-layer approximation (i.e., a zero normal 
pressure gradient). It is, however, a much milder approximation, since 
constant pressure is not applied through the boundary layer but over one 
grid line in the boundary layer. This approximation has yielded stable 
computation for both the non-separated and separated boundary layers 
[17]. For general curvilinear coordinates, Eq. (A. 25b) can be used to 
express Eq. (5.4) as 
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= 0 


(5.5a) 


This equation is approximated using a central-difference approximation 
on the wall surface, and a second-order backward-difference approxima- 
tion normal to the wall. 


P i,l,k * P i,2,k " P i,3,k + 2 ^ C 2 (P i+l,2,k ' P i-l,2,k ) 


+ C 1 (P i,2,k+l " P i,2,k-1 ) ^ /3 ’ 


(5.5b) 
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where 
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g 


(5.5c) 


gl2, g 23 anc j g22 are (j e fj ne( j by £q S> (A. 9)- (A. 10) . Then, density is 
computed based on pressure and wall temperature. At the wall, all 
boundary conditions are second order accurate and are satisfactory even 
for a skewed grid. Leading-edges of the wing have high curvatures near 
the back of the wing, therefore the normal pressure may not be equal to 
zero. Consequently, the above boundary condition may not be physically 
viable. However, the results indicate they are accurate enough for this 
problem. 

■ t ■ ;t t 

In the case of H-grid, a zero gradient in the q-direction is 
assumed for the symmetry boundary, i.e. 


au I 

5ti li.l.k 


= 0 . 


(5.6a) 


A backward-difference approximation is used to approximate the Eq. 
(5.6a) which results in 




(5.6b) 


Also, the velocity components normal to this boundary is set equal to 
zero 


w = 0 . 


(5.6c) 
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5.5 Boundary Conditions for Blunt Leading-Edge 


There are four boundaries in the computational domain with five 
different boundary conditions. They are upstream, downstream, outer, 
solid and symmetry boundary conditions. The top boundary (j=JL) 
contains the upstream and the outer boundaries. The upstream boundary 
condition is assumed to be the same as the freestream condition which 
can be expressed as 

V k ■ " u . • (5 - 7) 

Similar to the previous case, the outer boundary is located far away 
from the body to avoid any influence on the interaction region. Using a 
backward difference approximation, boundary condition for the top 
boundary can be expressed as 



= U 


JL-l,k * 


(5.8) 


A zero gradient in the C-direction is assumed for the downstream 
boundary. Using a backward-difference approximation, the following can 
be written 


U 


j.KL 



KL-1 * 


(5.9) 


At the wall, all velocity components are assumed to be zero. A zero 
normal pressure gradient is assumed for the solid surface, which can be 
expressed as 


'l.k ~ [4 P 2,k " P 3,k + 2 C 3 


(P 2,k+l ' P 2,k-1 )]/3 


(5.10) 
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where 


Co = 


T) C + T) C 

y y z z 
T7“2 — 

*1, + T)„ 


A zero gradient in the y-direction is assumed for the symmetry 
boundary, application of a backward-difference approximation yields the 
fol lowi ng 


U 


j.l 



(5.11a) 


Also, the velocity component normal to this boundary is set equal to 
zero 


v = 0 . 


(5.11b) 
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Chapter 6 


ADAPTIVE GRID GENERATION 

For hypersonic flow about blunt bodies (Fig. 2.3), the temperature, 
pressure and density of the flow Increases almost explosively across 

shock wave. At the same time, curved shock waves are close to the 
body. Numerical simulations of this phenomena have been a great 
challenge to computational fluid dynamics researchers. Presently, there 
is a great deal of interest in improving the quality of numerical 

\ * 

simulation techniques and grid adaption is one way to achieve this goal.. 

As stated earlier, grid generation is the first step in the 
numerical solutions of partial differential equations for complex 
geometric domains. Basically, grid generation is the creation of 
boundary-fitted curvilinear coordinates. The second step is the 

construction of difference equations for the partial differential 

equations. It is apparent that the accuracy of finite difference 
solutions depends on the fineness of the grid. Therefore, the finer the 
grid, the more accurate the numerical solution will be. Also, the 

accuracy of solutions depends on the resolution of the solution 
gradient. The presence of large gradients causes the error to be large 
In the difference approximation of derivatives. In the presence of a 
shock wave, more artificial diffusion must be added to retain adequate 
smoothness of the solutions. Therefore, there is a need for schemes that 
can resolve large gradients without adding additional grid points. An 
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adaptive scheme moves the grid points to regions of high gradients, when 
the locations of these gradients are not known a priori. Also, an 

adaptive method reduces the total number of grid points required to 
achieve a given accuracy, but it requires more computer time. In some 
instances, the computer time makes this method impractical. The ideas 
used in the construction of adaptive grid techniques are limited only by 
one's imagination; and any scheme that works in the sense of providing a 
better solution is a good one. The ultimate answer to numerical 
solutions of partial differential equations may well be to dynamically 
adapt grids, rather than to devise more elaborate difference 

representations and solution methods [11, 18], 

6.1 Literature Survey 

Adaptive methods have been used in the solutions of ordinary 
differential equations. In order to control the local truncation error 
[19], variable-step initial-value problems are solved by adjusting the 
step size as the integration advances. Adaptive methods have also been 
implemented for solving equations of motion in conjunction with the 
method of lines [20]. In this case, the time step is automatically 
adjusted to control local error. Similarly, adaptive methods have been 
used to solve boundary value problems [21-26]. An optimal grid for a 

two-point boundary value problem can be determined either implicitly or 
explicitly. In the implicit approach, the weight function depends upon 
the solution. As a result the original boundary value problem is 
converted into an augmented system in which the dependent variables and 
the grid are computed simultaneously. In the explicit approach, the 
weight function does not depend on the solution. Instead, it depends 
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upon a previously calculated solution. In the implicit approach, a 
nonlinear two-point boundary value problem should be solved even for a 
linear problem. Implicit techniques do not preserve the linear/non- 
linear character of the original problem. Moreover, even for the 
nonlinear problem, the augmented system is usually more difficult to 
solve than the original problem. On the other hand, the explicit 
technique preserves the linear/nonlinear character of the original two- 
point boundary value problem. 

Adaptive schemes are divided into two basic categories: 
differential and algebraic. Differential methods are based on a 
variational approach. Brackbill and Saltzman [27-30] have developed a 
technique for constructing adaptive grids using a variational approach. 
In their scheme, a function which contains a measure of grid smoothness, 
orthogonality and volume variation is minimized by using a variational 
principle. The smoothest grid can be generated by solutions of Laplace 
equations which are better known as elliptic systems [31-32]. This 
approach ignores the effects of orthogonality, and it is very slow. 
This method has been modified for better efficiency by dropping the 
second derivative terms in one coordinate direction [33]. This makes 
the equations parabolic, therefore they can be solved by a marching 
technique. A method which considers the orthogonality and volume 
variation has been developed by Steger and Sorenson [34]. This method 
is widely known as the hyperbolic method, and it can be solved by a non- 
iterative marching technique. The variational approach provides a solid 
mathematical basis for the adaptive methods, but the Euler-Lagrange 
equations must be solved in addition to the original governing 
differential equations. On the other hand, an algebraic method requires 
much less computational effort, but the grid may not be smooth. 
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Ra 1 and Anderson [35-40] have developed an algebraic technique 
where the grid movement is governed by estimates of the local error in 
the numerical solution. This is achieved by requiring the points in the 
large error regions to attract other points and points in the low error 
region to repel other points. Nakahashi and Deiwert [41] have formu- 
lated an algebraic method which is based on the variational principle. 
To reduce the overall solution error, a spring analogy is used to 
redistribute the grid points. In this case, operator-splitting and one- 
sided controls for the orthogonality and smoothness are used to make the 
method practical, robust and efficient. Ifc/yer [42-45] has used an 
adaptive method in which the points are moved along one set of the 
original coordinate lines in response to the evolving gradients in the 
physical solution. The analysis shows that the percentage change in a 
dependent variable can be determined a priori. Improvement in speed by 
an order of magnitude is obtained, but some problems with excessive 
skewness are encountered. 

Generally, dynamic adaption can be performed in two ways. One is 
to keep the computational space fixed and include the grid speed in the 
flow field equations. This is an ideal method to use for unsteady flow. 
The second way is to set the grid speed equal to zero and interpolate 
the solution onto the new grid after each adaption. The first way is an 
ideal method to use in unsteady flow while the second way is equivalent 
to solving a sequence of boundary-value problems and is an economical 
way to treat steady flows where solutions are approached asymptotically. 
In the second approach it is generally sufficient to adapt just a few 
times during the course of the computation. In this approach, the grid 
distribution at time N+l is determined from time N. Dwyer adapts the 
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grid points after each integration step or after a selected number of 
steps [42], However, the grid speed can be obtained by postulating a 
law governing it which is based on some solution properties. These 

equations can be integrated with the governing partial differential 
equations to yield the new grid distribution [46]. The advantage of 

this technique is that the location of grid points and grid speed are 

time accurate. 

The literature survey indicates that most techniques adapt to just 
one variable. This means that the weight function is based on the 

solution of one variable only. However, the solutions of equations of 
motion produce several dependent variables. Viscous-hypersonic flow 
over a blunt-body has large gradients in pressure, velocity, etc. in 
different parts of the flow field. For instance, there is a large 
gradient in pressure near the shock region, and at the same time there 
is a large gradient in velocity near the solid body. Therefore, there 
is a need for the development of an efficient grid adaption method which 
utilizes several variables simultaneously. 

6.2 Methods of Grid Adaption 

One reason to use grid adaption is to minimize the error over some 
domain by rearranging the grid. Calculus of variations can be used to 
perform this minimization. In general, a weighted integral, which is a 
measure of some grid or solution property over some domain, can be 
expressed as 

I = / W d V, (6.1) 

V 
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where W Is the weight function to be minimized. The selection of W may 
vary from problem to problem. There is a collection of definitions for 
W in [1], The weight function can be based on grid properties such as 
cell-volume, the average of the square of diagonal lengths, the cell 

area/volume ratio and cell skewness [33], There exists a differential 
equation which minimizes the integral I in Eq. (6.1). This differential 
equation is called the Euler-Lagrange equation [47]. The Euler-Lagrange 
equation can be found in [27-30], 

6.3 Multidimensional Grid Adaption 

Brackbill and Saltzman [27-30] have developed a technique based on 
a variational approach. In their scheme, a function which contains 
measures of grid smoothness, orthogonality, and volume variation is 

minimized. To maximize the smoothness of the grid, the following 

integral must be minimized 

3 

I = / z V S 1 • V S 1 dV . (6.2) 

s y i=l 

This is simply the sum of the squares of cell-edge lengths. Similarly, 
orthogonality can be acquired by minimizing the integral I 0 

I » / (7 C 1 • 7 f*) 2 dV , (6.3) 

o v 

with (i,j,k) cyclic. This integral vanishes for an orthogonal grid. 
The concentration or cell-volume variation can be obtained by minimizing 
the integral 


I w - I W 0 dV , 
w v 


(6.4) 


43 



where W is a specified weight function. This causes the cells to be 
small where the weight function is large. The grid generation system 
which provides smoothness, orthogonality and concentration is obtained 
by minimizing the total integral I which is a linear combination of I s , 
I 0 and I„ 

1 ’ h + >o * \ \ ■ (6 ' 5) 

The competing features such as smoothness, orthogonality, and cell 
volume variation can be stressed by the proper choice of the coeffi- 
cients X Q and \ w . For example, a large \ Q will result in a nearly 
orthogonal grid at the cost of the smoothness and the concentration. 
The Euler-Lagrange equations for the sums of those individual integrals 
form the system of partial differential equations from which the 
coordinate system is generated. The equations are quasi-linear, second- 
order partial differential equations with coefficients which are 
quadratic functions of the first derivatives [29]. This variational 
formulation is equivalent to Winslow's method [31] where X and X,, are 
set equal to zero. The Euler equations are those given by Winslow, and 
their solution maximizes the smoothness. This is also used by Thompson 
et al. [32]. The additional terms alter other characteristics of the 
mapping in a similar way. The cell-size variation and skewness can be 
controlled by proper selection of I $ , I Q and I w . The use of a 
variational approach provides a solid mathematical basis for grid 
adaption. But, the Euler-Lagrange equations must' be solved in addition 
to the governing equations of fluid motion. For further information, 
readers are refered to excellent articles by Thompson [11, 18]. 
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6.4 One-Dimensional Grid Adaption 


The Euler-Laorange equation Eq. (6.5), is general and capable of 
adapting grids simultaneously in multiple dimensions. When the solution 
varies predominately in a single direction, one-dimensional adaption can 
be applied with the grid points constrained to move along one family of 
fixed curvilinear coordinate lines. The fixed family of lines is 
established by generating a full multidimensional grid using any 
standard grid generation technique. The points generated for the 
initial grid together with some interpolation procedure, e.g., cubic or 
linear interpolation, serve to define the fixed lines along which the 
grid points will move during the adaption. This is done explicitly, 
therefore there is no need to solve any differential equations. 

A technique called equidistribution is developed to improve the 
solutions of boundary value problems [21-26]. This technique has proven 
to be effective and efficient. This technique is used to minimize the 
error by redistributing grid points such that a weight function is 
constant over each interval. The Euler-Lagrange equation is 

x^ W = constant . (6.6) 

This minimizes the following integral 

I. = / W(£) xf d£ . (6.7) 

1 0 * 

Equation (6.7) represents the energy of a system of springs with the 
spring constant W(5). spanning each grid Interval. The weight 
function is associated with the grid points themselves and not with 
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their locations. An alternative viewpoint results from integrating over 
x. Instead of £, i.e., summing over the grid intervals rather than over 
the grid points. This can be expressed as 

1 c 2 

h • ! 0 -HTXT dx • (6 - 8 > 

The Euler-Lagrange equation for this formulation is given by Eq. (6.6). 
£ x is considered to represent the point density. This variational 
problem represents a minimization over the density of the grid points 
subjected to a weight function. This can produce smooth grid distribu- 
tions. Here the weight function W(x) is associated with the location of 
grid points. If the weight function is associated with the grid points 
themselves rather than their locations, W = W(£). Equation (6.6) is 
the Euler equation for the following integral 


. f 1 r ] 2 

j 0 L *irr J 


d£ 


(6.9) 


where l, is a measure of the smoothness of the grid distribution, with 
the emphasis placed on smoothness in certain regions. This is inversely 
proportional to the weight function W(£). Equation (6.6) is the Euler- 
Lagrange equation for the integrals in Eqs. (6. 7-6. 9) which can be 
written as 


x 

5(x) = / W(t) dt . (6.10) 

0 


Equation (6.10) can be written in terms of arc length as 


s 

$(s) = / W(t) dt . (6.11) 

0 
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The weight function is used to reduce the grid point spacing where W is 
large, and the weight function should be some measure of the error. 
White [26] has suggested the following form of the weight function 

H = [a + |U (m) | 2 ] , (6.12) 

where a is constant. With m=l and a=0, this becomes 

W = |U | . (6.13) 

A combination of this equation and Eq. (6.6) yields, 

= Constant . (6.14) 

This choice replaces the grid points so that the same change in the 
solution occurs at each grid interval. This is simply the solution 

gradient. Taking n=l and o=l yields 

W = r i + |uj 2 . (6.15a) 

Combination of Eq. (6.15a) and Eq. (6.6) results in 

J~~7 T 

X C + = S C = Constant • (6.15b) 

This produces a uniform distribution of arc lengths on the solution 

curve. White's results [25] indicated that the arc length form is 
favored. The disadvantage of this method is that the weight function 
near the extreme solution, i.e. U x =0 locally, is treated as a flat 

region. Concentration near the solution extreme can be achieved by 
incorporating some effect of the second derivative (U xx ) into the 
weight function [42] such as 
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(6.16a) 


W = 1 + a f(U x ) + p g(U xx ) , 

where a and p are positive parameters. Equation (6.10) can be 
rewritten in normalized form as 

x 

/ W ( t ) dt 

5(x) --j . (6.16b) 

/ W(t) dt 
0 

With the second derivative terms included, the value of a must be 
continually updated to keep the same relational emphasis or concentra- 
tion. Therefore, a system of two equations and two unknowns must be 
solved for each fixed grid line [43-45]. It will be shown later that 
through the reformulation of Eq. (6.16b), the parameter can be found 
directly or the matrix inversion can be entirely avoided. 

6.5 Cne-Dimensional Grid Adaption With Several Variables 

Flow equation solutions consist of several variables. Therefore, 
the weight function should also be a function of more than one variable. 
It is desirable to devise a scheme in which grid points can adapt to 
several variables with control of the magnitude of adaption for each 
variable. In the case of high-speed flow, velocity has large gradients 
in some regions where pressure is constant or vice versa. In general, 
the weight function can be expressed as 

N 

W = 1 + E b. f. (6.17) 

1=1 1 1 

where N is the number of variables, b^ are constants, fj are variables 
or their derivatives, and 1 is for uniformity. .A substitution of Eq. 
(6.17) into Eq. (6.16b) results in 


i 
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5(S) = 


f 


(6.18) 


N 

S + £ b. F, (S) 

1 = 1 1 1 

S max + b i F i * S max^ 


where 


F. (S) 



(t) dt . 


It should be noted that Eq. (6.18) is for adapting along a fixed grid 
line. To ensure that £(S) increases monotonically, and f^ should be 
positive. In order to keep the same relative emphasis of the concen- 
tration along each grid line, b^ should be computed based on some 

percentage of the grid points being allocated to each variable. The 

percentage of grid points assigned to a particular function f n - can be 

expressed as 


R. = 

J 


b j F j ^ S max* 

— R 


, j=l,2,...,N . 


(6.19) 


^max + b i F 1 ^max^ 


Rearranging this equation results in 


C V { V ' {c j> 


where 



F. (S ) 
D l max 

= R j ITT S— 7 


J 


max 


1 * J 


( 6 . 20 ) 


R. - 1 
J 


1 = ■ - R j W F j <W 


Therefore, a system of N equations and N unknowns must be solved for 
each fixed grid line (N is the number of variables). This can be 
avoided by the reformulation of Eq. (6.18). The crucial steps are 
outlined here. Equation (6.19) can be rewritten as 
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( 6 . 21 ) 


b, = r — — r [S + Z b. F. (S )] 

J F 7 T! w T max 1-1 1 1 max 


A substitution of Eq. (6.21) Into Eq. (6.18) results In 


C(S) = 


N R F. ($) 

S + W max . E . TTTS J 

i=l i max 

FT 

$max + ''max 


where 


W = [s + Z b. F. (S )1 
max L max ._j i i ' max'- 1 


A summation of Eq. (6.19) over all j values yields 


N 

m Z b. F. (S ) 

" i_i J J max 7 

Z R. = — 

J 


J-l 


T 


s + Z b. F. (S ) 
max , =1 i i v max' 


Rearrangement of Eq. (6.23) results in 


N 

Z 

1=1 


F. (S ) 
i max 


N 

S Z R. 
max i 

R 

1 - Z R. 
i=l 1 


A substitution of Eq. (6.24) into Eq. (6.22) yields 


5(S) = 


[1 - 


max 


N 

z 

1=1 




N 

Z 

1=1 


R. 

l 


F i (S) 

=773 7 

1 max 


( 6 . 22 ) 


(6.23) 


(6.24) 


(6.25) 
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This reformulation avoids the need for continuous updating the b-j's to 
keep the same relative emphasis on concentration. Applying this 
equation, grid points can be adapted with more than one variable without 
any need for matrix inversions. Eiseman [48] proposed a very similar 
approach but did not elaborate on this. 

Presently, Eq. (6.18) is approximated by a trapezoidal rule. For 
unequally-spaced data, this can be written as 


S . JL 

F. (S) = / f. (t) dt =4 E C(t. , + t.) (f.(t.) 
l 0 i ^ j=l 3 1 3 1 J 



D) , 


where 

F.(0) = 0 . 

In the initial stages of the solutions, there exist large oscilla- 
tions in the flow and fluid properties. Consequently, the adapted grid 
will have these oscillations as well. In order to have a smooth grid, 
these oscillations can be smoothed out with the following filter 



/r n + + f n + f n )/5 

r jk + r j+l,k + r j-l,k + r j,k-l + r j,k+l' /b 


where 

r = (y,z) T . 


(6.25a) 


This is equivalent to the Lapalace filter which can be expressed as 

5 2 r + afr = o # (6.25b) 


5ti 5 C 

This slightly reduces the effect of adaption, but it filters out low 
frequency oscillations. 
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In some physical phenomena like shocks, the length of discontinuity 
is of the order of molecule size. If the weight function contains the 
derivatives with respect to physical coordinates, the adapted grid tends 
to have a very small spacing. This can be corrected by reducing the 
effect of the weight functions near the discontinuities. Therefore, the 
weight functions (f^) can be multiplied by the following function 

AS 

' AT~ 

(1 - e min ) 

where AS . is some allowable minimum spacing. This function varies 
from zero to one and is proportional to the spacing. Also, this can be 
corrected by replacing the derivatives of physical coordinates with the 
derivatives of the computational coordinates. 

After the new adapted grid has been created, the solutions are 
interpolated into the new grid distribution. Presently, a piece-wise 
linear spline has been used. However, this may create some problems in 
the case of unsteady flow. Either higher order interoplation should be 
used or the time derivatives of the grid point should be incorporated 
into transformed governing equations. 
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Chapter 7 


RESULTS AND DISCUSSION 

In previous experiments [9], the surface oil-flow patterns over the 
Butler wing at various angles of attack and at Mach number 3.5 show no 
signs of transition and the nature of the oil streak lines is typical of 
a laminar flow. Therefore, results are obtained only for laminar flow 
over a Butler wing at a Mach number of 3.5, Reynolds number of 2xl0 6 
(based on chord length), free stream static temperature of 390°R 
(216. 67°K) , wall temperature of 1092°R (606.67°K), length of 0.80 ft 
(0.2438 m), specific heat ratio — of 1.4, specific heat at a constant 
volume of 4290 (ft/sec)^/R, and at zero and ten degrees angle of attack. 
In this study, a two boundary grid generation (TBGG) technique [13] is 
used. This method is essentially an algebraic method. The application 
of the TBGG method requires that the entire body be sliced into 
different cross sections. These cross sections are obtained in the 
stream-wise direction by analytical descriptions of the wing surface, 
(Eqs. 2.1a-2.1c). Then, two types of grid are generated for this wing; 
the H-type and 0-type. Then, results of both cases are compared and 
discussed. 


7.1 H-Type Grid 

In this case, the entire flow field is sliced Into fifty-five 

stations in the stream-wise direction, and each station has 64x36 grid 

points (Fig. 7.1). There Is a total of 126,720 grid points which take 
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Fig. 7.1 H-type grid for a Butler wing 
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Fig. 7.1 Continued 
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Fig. 7.1 Continued 
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Fig. 7.1 Continued 



2.8 million 32-bit words of primary memory (16 variables). The required 
computational time is 1.9xl0 -5 sec/grid-point/iteration (2.5 
sec/iteration) which is typical for a CYBER 205 with two pipes. Results 
are obtained for zero angle of attack. The computed pressures are 
plotted in (Figs. 7.2a-7.2c). The pressure coefficient along the center 
line is shown in Fig. 7.2a. The results are compared with available 
experimental and numerical results [6-9]. The results on the center 
line are in excellent agreement with the experimental and previously 
obtained numerical results (Refs. 6-9). At 41.66? and 68.33? chordwise 
position, the pressure ratios are plotted against the conical span-wise 
coordinates y/xtan(e), (Figs. 7.2b-7.2c). They are in good agreement 
with experimental and numerical results. However, there are some 
discrepancies in the results between 30° and 60°. This is probably due 
to the fact that grid lines are not orthogonal near those regions and 
this may be a direct consequence of the H-type grid. Also in this case, 
the wing's leading edges are represented by a jagged-line, in other 
words grid lines are not along the leading edges of the wing. This does 
not allow us to compute the boundary conditions using a second-order 
accurate formula. Therefore solutions are not second-order accurate 
near solid boundaries. 


7.2 0-Type Grid 

For this case, the physical domain is limited to 5? to 95? of the 
wing. This is done to avoid any singularities. The conical Navier- 
Stokes solutions are enforced for the upstream boundary which is located 
at 5? of the wing. This solution is obtained by the integration of the 
Navier-Stokes equations [16] for a conical grid with proper boundary 
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AXIAL DISTANCE ( X/C) 


Fig. 7.2a Pressure coefficient along the center line 



ANGLE (X/C = 41 .677.) 


Fig. 7.2b Pressure ratio at 41.67% of chord 



conditions. The wing is sliced Into forty-one stations in the stream- 
wise direction, and each station has 41x127 grid points (Fig. 7.3). 
There is a total of 213,487 grid points which take 4.7 million 32-bit 
words of primary memory. Results are obtained for zero and ten degrees 
angle of attack. 

Results for zero angle of attack are compared with results obtained 
with the 0-type grid, with results from the experiments and other 
numerical results [6-9], The computed pressure is plotted in Figs. 

7.4a-7.4c. The pressure coefficient along the center line (Fig. 7.4a) 
is in good agreement with other numerical and experimental results. 
Nevertheless, there is some discrepancy near the nose region. This may 
be due to the fact that the upstream solutions are based on the conical 
solutions. But, these solutions match exactly with results from H-type 
grids. This is because the grid topology near the center line is the 
same for both grid types. At 41.67* and 68.33* chordwise positions, the 
pressure ratio is plotted against the conical span-wise coordinates 
y/xtan(6). They are in excellent agreement with experimental and 

numerical results (Figs. 7.4b, 7.4c). In addition they are much closer 
to the experimental results compared to the results obtained from the H- 
type grid. This may be due to good grid orthogonality in the case of an 
0-type grid, and the wing's leading edges are represented by a straight 
line. On the thick sections near the nose the pressure is highest on 
the centerline and falls toward the leading edge. Figure 7.4d shows the 
cross-flow velocity at 5*. 23*, 41*, 59*, 77* and 95*. 

The results for ten degrees angle of attack are compared with 

experimental results. The computed pressure is plotted in Figs. 7.5a- 

7.5d. At 17*, 30*, 50* and 70* chordwise positions, the pressure ratio 
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Fig. 7.3 Continued 
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Fig. 7.3 Continued 
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Fig. 7.3 Continued 


PRESSURE RATIO 



Fig. 7.4a Pressure coefficient along the center line 



Fig. 7.4b Pressure ratio at 41.67% of chord 





Fig. 7.4d Cross-flow velocity vectors 
(zero angle of attack) 
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Fig. 7.4d Continued 


70 










Fig. 7.4e Streak lines (zero angle of attack) 
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Is plotted against the conical spanwise coordinate y/xtan(9). They are 
in good agreement with experimental and numerical results (Figs. 7.5a- 
7.5d). As in the previous case, on the thick sections near the nose the 
pressure is highest on the centerline and falls toward the leading edge, 
whereas near the trailing edge the spanwise distribution is more "wing 
like" with the maximum pressure at the leading edge. The change-over is 
shown by the pressure peaks in the pressure distributions at x/c=0.5 and 
0.7 at 10 degrees angle of attack. There are some discrepancies near 
x/c=30?-50?; this may be due to the fact that Squire [9] has not used 
the exact model of the Butler wing. In order to mount the model in the 
wind tunnel, the lower surface is distorted to include a sting 
support. Figure 7.5e shows the cross-flow velocity at 5?, 23?, 41?, 
59?, 77? and 95?. These figures show a weak cross-wise separation on 
the suction side which is confined to the body. At 59?, the cross flow 
has separated but a well-defined vortex is not visible. Squire [7] has 
performed a series of tests to investigate the effects of thickness on 
the longitudinal characteristics of a delta wing with different aspect 
ratios. The tests on the thick symmetrical delta wings have confirmed 
that the lift curve slope decreases as the thickness increases. This 
loss of lift is associated with a weaker vortex system giving less 
nonlinear lift. Squire [9] also observed a pair of vortices at the 
trailing edges, but there was no sign of any span-wise flow outboard of 
these vortices. The Butler wing has a round leading edge for most of 
its length, and previous numerical experiments have indicated that flow- 
field solutions are inconsistent using both the Euler and the Navier- 
Stokes equations for this type of geometry [49]. The problem appears to 
be the determination of the initial location of separation over the 
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Fig. 7.5a Pressure coefficient (x/c=17% 
ten degree angle of attack) 
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Fig. 7. 5d Pressure coefficient (x/c=70% 
ten degree angle of attack) 




Fig. 7 . be Cross-flow velocity vectors 
(ten degree angle of attack) 




Fig. 7.5e Continued 





Fig. 7.5e Continued 





Fig. 7.5f Streak lines (ten degree angle of attack) 
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Fig. 7 . b f Streak lines (ten degree angle of attack) 
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suction side. If the leading edge is sharp, there Is no problem since 
the separation line is along the leading edge. For a rounded leading- 
edge, the separation point is not defined uniquely because the numerical 
dissipation influences where the separation point occurs [49]. The 
leading-edge separation for the rounded-edge may be a numerical 
phenomenon. It Is possible that this problem of uniqueness occurs in 
this study as it does with all other numerical experiments, having 
rounded leading edge and artificial dissipation applied in the solution 
technique. 

The question has been raised concerning the application of the 
Navier-Stokes equations in a flow field study when only the pressure 
solution and streak lines are compared with experimental data. It is 
argued that the solution of the Euler equations would produce the same 
information with a high degree of accuracy. This is a legitimate 
argument if the sole purpose of the study Is to collect inviscid 
Information. However, we believe that the solution of the Navier-Stokes 
equations provides considerably more information if there are enough 
grid points in the proper locations. Also, the solution of the Navier- 
Stokes equations provides experience for the future when computer speed 
will be fast enough to handle the large number of grid points. 

7.3 Static Grid Adaption 

The second part of this study concentrates on finite-difference 
methods in which the grid points adapt to the solution dynamically to 
obtain an accurate solution for hypersonic flow. A computer program has 
been written to utilize Eq. (6.25) for grid adaption. Presently, this 
code is being run on the Network Operating System (NOS) and the Vector 
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Processing System (VPS-32) at NASA Langley Research Center. Hypersonic 
flow over a blunt nose is a typical test case in computational fluid 
dynamics. This problem has a detached shock which should be resolved 
accurately, and the location and the magnitude of the shock are not 
known a priori. The grid should adapt as the solution progresses. This 
problem is used to analyze and verify the adaptive method. The 
equations cf motion are solved by the MacCormack method [12-13] for 
hypersonic flow over a small-radius blunt-body with the inclined-plate 
afterbody (Fig. 2.3). The blunt leading edge is a part of the panel 
holder which has been tested at the Langley Research Center [50], The 
results are obtained at the following conditions: free stream Mach 
number of 6.8, a pressure of 9.26 lb/ft 2 , velocity of 6510 ft/sec, 
temperature of 375°R (static temperature), Reynolds number of 220,000, 
specific heat ratio of 1.38, universal gas constant of 1771 ft 2 /sec 2 /R 
and a wall temperature of 540° Rankine. 

Two tests have been performed: static adaption and dynamic 

adaption. For static adaption, the solution has been obtained with the 
fixed grid points shown in Fig. (7.6). Then, the grid points are 
adapted to two variables, the first and the second derivatives of 
pressure. Results are shown in Fig. (7.7). In this case, twenty 
percent of the grid points are allocated to first and second derivatives 
of pressure (R 1 =R 2 =202). For the same case. Fig. (7.8) shows adaption 
with R2=R2=50$. In this case, all the grid points are allocated to the 
first and the second derivatives of the pressure. This explains the 
large voids in the constant pressure regions. Figures 7. 7-7. 8 lack grid 
resolution in the vicinity of the solid boundaries. This is due to the 
constant pressure near the solid boundaries. But Eq. (6.25) can adapt 
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Fig. 7.6 Initial grid distribution 
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Fig. 7.7 Adapted grid (R^=R2=20%) 
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Fig. 7.8 Adapted grid (R-j=R 2 =50%) 
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to several variables. Figure 7.9 shows the grid points which are 
adapted to two variables — pressure and velocity. The weight function 
consists of the first derivatives of pressure, the velocity and the 
second derivative of pressure. Twenty percent of the grid points are 
allocated to each function, and forty percent of the grid points are 
allocated for uniformity. This avoids the creation of any void. It is 
noted that grid points are clustered near the shock and the solid body. 
Therefore, it is possible to resolve the pressure as well as the 
velocity gradients in the boundary layer region. In Fig. (7.10), thirty 
percent of grid points are allocated to the first derivative of Mach 
number. If there was a chemical reaction in process, some of the grid 
points could have been allocated for resolving the gradients of the 
chemical species. 


7.4 Dynamic Grid Adaption 

The above procedure has been applied dynamically to the same 
problem. Figure 15a shows the initial grid distribution for this 
problem. Figure 7.11 shows sequences of grid distributions at a 
different time. In this case, grid points are adapted to six variables: 
the first and second derivatives of pressure, Mach number and 
velocity. Ten percent of the grid points are allocated equally to first 
derivatives, and five percent of the grid points are allocated to their 
second derivatives. Fifty-five percent of the grid points are also 
allocated to the uniformity. A movie has been produced of this work 
which shows the dynamic adaption. A few frames are shown in Fig. 
(7.12). They demonstrate how grid points are attracted toward high 
gradient regions and repelled from low gradient regions. 


3 / 
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Fig. 7.9 Adapted grid (R-j =R 2 = R 3 = 20 %) 
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Fig. 7.10 


Adapted grid (Rg=30%) 
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Fig. 7.11 Adapted grid (Dynamic) 
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Fig. 7.11 Continued 
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Fig. 7.11 Continued 
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Fig. 7.11 Continued 
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Fig. 7.12 Adapted grid (dynamic) 
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Continued 
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Fig. 7.12 Continued 
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Chapter 8 


CONCLUDING REMARKS 

General formulations are presented to investigate the flow field 
over complex configurations for high-speed free stream conditions. An 
advanced algebraic method is used to generate grids around these 
configurations. The computational procedure developed is applied to 
investigate the flow field over a Butler wing. Illustrative results 
obtained for specified free stream conditions compare very well with 
available experimental and numerical results. Results are obtained for 
laminar flow over a Butler wing at a Mach number of 3.5, Reynolds number 
of 2xl0^/ft (6.56xlC^/m) , free static temperature of 390°R (216.67 K), 
wall temperature of 1092°R ( 606 . 6 7 ° K ) , length of 0.80 ft (0.2438 m) when 
the wing is at zero and ten degrees angle of attack. Two types of grids 
have been generated for this wing; H-type and 0-type. Results of both 
cases are compared and discussed. The future plans to extend this study 
are to include the sting, use multiple grid, different Mach numbers, and 
to include turbulence. 

A grid adaption method has been developed with the capability of 
adapting grid points to several variables. This method is an algebraic 
method, and has been formulated in such a way that there is no need for 
any matrix inversion. The method is used in conjunction with the 
calculation of hypersonic flow over a blunt-nose body. A movie has been 
produced which shows simultaneously the transient behavior of the 
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solution and the grid adaption. The results indicate the viability and 
validity of the proposed method. The future plans for this technique 
are to use this problem with a true unsteady problem, to study the 
effect of interpolation, to consider more complex geometries and finally 
to adapt the problem in three-dimensions. 
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APPENDICES 


APPENDIX A 


MATHEMATICAL DETAILS FOR THE TRANSFORMED EQUATIONS 
A.l Curvilinear Coordinates 

In covariant coordinates system (x.j)» the position vector of a 
point from the origin is expressed as 

r = e. x. - e 2 xj + e £ x 2 + e 3 x 3 (A.l) 

where e^ is the covariant base vector. 

In the present study, covariant coordinates are labeled as x, y, 
and z (i, 3, k) and contravariant coordinates are labeled as S, n and 
C (i, J, k). The covariant base vectors are defined as 



where J is the Jacobian of transformation. Magnitude of Jacobian (|j|) 
is the local value of the ratio of an elemental volume in the mapped 
(usually cube) cell to the corresponding elemental volume in the 
physical (usually distorted) cell. 
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The contravarfant base vectors are defined as 



Position vector can be expressed explicitly in 
contravariant vector (x^); however, the infinitesimal vector 
expressed as 


dr = XX d x* = e. dx 1 ' . 
dx 

Also the magnitude of arclength (ds) can be expressed as 


( ds ) ^ = dr • dr 


Substitution of Eq. A. 4 into Eq. A. 5 will result in 


ds z = (e. • e.) dx . dx . = g. . dx. dx . 

I J ' J • J * J 


where g^j is called covariant fundamental metric coefficients. 
These coefficients can be defined as 


_ n -l-,T n _1 n d k 5 k 
J dx dx J 


They are defined as 


'll 


2 2 2 

h + h + z f 


g l2 = 921 = x n + y S y n + z c z ri ’ 


(A. 3) 

terms of 
dr_ can be 

(A. 4) 
(A. 5) 
(A. 6) 

(A. 7) 

(A. 8a) 
(A. 8b) 


117 



9 


(A. 8c) 


g 13 = 9 31 *1 X C + y C + Z l Z i 


2 2 2 
g 22 = X T! + y r, + z n 


g 23 = g 32 X n X C + Y n y C + Z ti Z C * 


2 2 2 
g 33 = X C + y C + Z C 


(A.8d) 

(A.8e) 

(A.8f) 


Similarly, contravariant fundamental metric coefficients are defined as 

(A. 9) 


g ij = e 1 • e j = [J] [J] T * 9x1 - 9xJ 




or 


11 2 2 2 

> -4 + Zy + 4 ’ 


21 _ 12 


g ^x + S' S + S ^z ' 


13 _ 31 


22 


g * «x Sc + S S + S S ’ 

n 2 + T) 2 + r| 2 , 

x y z 


23 .32 


.33 


g = \ C x + S Sr + ^ C z * 

c 2 + c 2 + c 2 . 

x y v z 


(A. 10a) 
(A. 10b) 
(A. 10c) 
(A.IOd) 
(A.IOe) 
(A.IOf ) 


Furthermore, there exists a unique relationship between contra- 
variant and covariant fundamental metric coefficients. 


ij g rs g *s ” g rt g ts _ 


fg 


ij 


where 


G 11 g 22 g 33 


G * . 

Tgfri 

,y ij ! 


2 

g 23 


Cofactor of g.. 

* J 




(A. 11) 


(A. 12a) 


G 12 = G 21 = g 13 g 23 " g 12 g 33 * 


(A. 12b) 
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G 13 = G 31 “ g 12 g 23 " 9 13 g 22 


9 


(A. 12c) 


where 


G 22 9 11 g 33 " 9 13 • 

G 23 = G 32 = 9 12 9 13 ~ g 23 9 11 ’ 
2 

G 33 = 9 11 9 22 " 9 12 ’ 




Cofactor of 




(A. 12d ) 
( A. 12e ) 
(A. 12f ) 


There Is also a relationship between covariant and contravariant base 
vector 


where 


e = 


£j x ik 


hjl 


17? * 


9i.il = l J-1 | 2 • 


(A. 13) 




i .e. 


or 


i ( Vc' y eV 1 ' ( Vc‘ x cV J + lx -i y <;' x cV k 

e ' ~FT, 


J 



J 


1 

FTi 


(y ^ z c- y c z T,> 

( Vc' y cV 


(X Z--X-Z ) 
t) C C n 


-(x.z-xz) 

£ T] T) c 


l Vc'Vi ) 


(x^z^-x^z^) -(x r y r -x r y P ) 


E, J C> ^ 


(x 5V\V 


where 


J = e 
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There Is also a relationship between contravariant and covariant base 


vector 


where 


i .e. 



x e k 


ii = 


9 


1J 


^T 



|g lj l - M 2 . 


(A. 15) 

^y^z^z^y) 

-KWx 1 

(rix^y ) 


- (5 yVW 

< 5 x^- 5 z C x ) 

- ( WyW’ 

(A. 16a) 

^y^z'^z^y ^ 

■ ({ xV^’ 




or from Eq. (A. 2) 


[J" 1 ] = Ce, 


*/■ 


w 


(y^zV '•tyWy 1 "yVW 

( Vx-Vx> ' (C xV 5 zV ? d J l • <A - 16bl 

VyV * 1 -‘'xVVx 1 U xVV>’ 


The relationship between vector bases can be obtained also by matrix 
algebra. From basic matrix identity, Jacobian can be written as 

T 

(A. 17) 


[J] = [J -1 ] = 


Transpose of cofactor [J *] [[J *] ] 


|J 


FT 


J 


FT 
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Equation (A. 17) is the same as Eq. (A. 14). There also exists an inverse 
relation for Eq. (A. 17). 


A. 2 Vector Representation in Curvilinear Coordinates 

A vector F(F x i, F y j, F z k) can be expressed in contravarient 
coordinates as 


F i ■ 


(A. 18) 


This is a projection of F^ on x^ coordinate 


where 


i 


i ii r i _ ax 1 r , ax 1 f . ax 1 r 
l e I F - F X + w F y + W F Z 


(A. 19) 


Equation (A. 19) can be expressed as 


V (s n> 


F^/(9n ) 


J 22 

J 33‘ 




[J] 




(A. 20) 


where is defined in Eq. (A. 8). 

The Inverse relation to Eq. (A. 20) is 


( V(si! /2 ) 

H j y(»* *) 



(A. 21) 
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For example, the velocity vector In covariant coordinates can be written 


as 



f U/(gj} /2 ) 

W 1 < v/(g 2 | /2 ) 

l W/(g 3 l /2 ) 


(A. 22a) 


Also, the velocity vector in contravarlant coordinates can be expressed 
as 


lu/(g n ) 
\ v/(g 22 ) 

( W/ ^ g 33^ 




(A. 22b) 


where u and U are velocities in the covariant and the contravariant 
coordinate systems, respectively. 

A. 3 Normal Derivative in Curvilinear Coordinates 
A normal derivative of a scalar variable (A) can be computed as 

where 

VS e i e i 

n = T^T = TeTT ""(TV 7 ^ ’ 

VA = A x i + A y j + A 2 k 


(A. 23a) 


(A. 23b) 


(A. 23c) 
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or 



Equation (A. 23a) can be written as 


9A 


constant 4> 



For a constant l, Eq. (A. 24b) can be written as 
8A | 9 U A 5 * g 12 A^ + g 13 A ? 

™ h ( 7 1 ) 172 


For a constant ti, Eq. (A. 24b) can be written as 


8A | S 21 \ * S ?2 ^ ♦ 9 23 A c 

(T 7 ? 72 


For a constant C, Eq. (A. 24b) can be written as 
aA | 9 31 a 5 ♦ g 32 A„ ♦ g 33 A c 

m k (T 7 ? 72 


(A. 24a) 


(A. 24b) 


(A. 25a) 


(A. 25b) 


(A. 25c) 
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also 


(A)' = 


i e 


V A = - 4 - Z (e 1 • eh A? 


e 1 ! |e | j=l 


i ; jj »a_ 

I^A 9 


(A. 26) 


e.g. 


f A >* ■ jTTp77 I 8 ” A C + g12 \ * 3 ' 3 \ ] • 


(A. 27) 


This equation Is similar to Eq. (A. 25a) 


A. 4 Miscellaneous Relations 


The angle between two grid lines is given by 


e 1 * e j 

cos e ij = KT1?7’ 


COS 9. . = 


i J 


,J tM | 9jj |]^‘ 


(A. 28) 


Therefore, for the orthogonal grid, the following should be true. 


=0 for i * j . 


(A. 29) 


Arclength is defined as 


(ds) = 


33 i j 

E E g . . dx dx J . 

1=1 j=l J 


(A. 30a) 


An arclength along x 1 coordinate is defined as 


(ds) 1 = (g ii ) 1/2 dx 1 . 


(A. 30b) 
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The area of an element on which the x 1 constant is defined as 


(dr 1 ' = |e. d • e. d x k | = (g..) dx 1 ' dx^ , 

J K I J 


e,g, 


y /gTT 

dl* = |e 2 dn x e 3 dc| = /EJJ dn dC = — ip dn dC , 


dr 11 = G 22 dl dC = — d£ dC , 


i 22 


dr c = d£ dn = — -p d£ dn . 


J 33 


The volume of an element is defined as 


dV I 0 ' 1 ! d5d " « 


= | e 2 * (e 2 xe 3 )| dl dn dC = (|g..|) 1/2 dl dn dC 


iJ 1 


(A. 31) 


(A. 32a) 


(A. 32b) 


(A. 33) 


125 



APPENDIX B 


TIME-STEP ANALYSIS 

In order to analyze the governing equations, they can be divided 
into two parts: Inviscid and viscous. This makes the task of analyzing 
much simpler. However, this is another assumption which may impair the 
resul ts. 


B.l Inviscid Part 

For the inviscid part, the governing equations can be written as 


P t + v • (pu) = 0 , 

(B.la ) 

(pu) t + V • (puu) + (V • P) 6^ . = 0 , 

(B.lb) 

e t + V • (ev) + [v • (pu)] 6^ = 0 , 

(B.lc) 

? ? ? 

P . p(u + V + w ) 
e -ipr + 2 ’ 

(B. Id) 

y P = pC 2 . 

(B.le) 


After some algebraic manipulations, Eqs. (B.la)-(B.ld) can be written as 

q t + A q x + B q y + C q z = 0 (B.2) 

where 


q 


(p, u, v, w, 
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Equation (B.2) can be transformed from physical coordinates to computa- 
tional coordinates as 


Q t + A + B + C = 0 , 

where 

A = A £ x + B + C l z , 

B = A ri x + B riy + C n 2 , (B.3) 

C - A C + B C + C C . 

a y i 

This equation can be split as follows 

," +1 -, n *ut 5 fl" E =o, 

Q n+2 . q n*l + B 4 t q q" +1 = 0 , (B.4) 

n+3 n+2 j. . . n+2 n 

q -q +CAt^q^ =0. 

This splitting is valid and stable if the time-step of each operator 
does not exceed the allowable step-size for each operator. It is 
consistent if the sum of the time steps for each operator be equal. 

This method would be second order accurate if the sequences of the 
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operator are symmetric. In order to analyze this equation, finite- 
Fourier series can be introduced to monitor the growth and decay of the 
error as 

q = qe st e m . (B.5) 

Substituting this equation in Eq. (B.4) results in 


e SM - 1 + A * 

A t. * 


ikAq _ -ikAn 


= 0 


(B.6) 


2 A ti 


Rearranging and collecting terms, yields the following 
A . At, 

e sAt = [G] = [I - i sin e)] , 0 = k A l . (B.7) 

The system is stable if the largest eigenvalue is less than unity. This 
condition Insures that the error always decays. 

| G | < 0 . (B.8) 


Therefore, the determinant of [G] can be set equal to zero 

_ At, 

det [a sin 9 - XI ] = 0 , 


or 



Pi E 

P5yE 


0 

Eu-X 

0 

0 

E5 x/p 

0 

Eu-X 

0 

E5 y/P 

0 

0 

Eu-X 

E5 ./p 

P« X E 

P« y E 

pcC z E 

Eu-X 
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where 


At- 

E = sine (A£ x + B C y + C C z ) » 


u = u l x + v l y + w l z . 


The determinant can be written as 


( u E - X ) 3 [(u E - X ) 2 - E 2 c 2 U 2 + Sy + l\)} = 0 . (B.10) 


Equation (B.10) can be solved for X as 
At. 


1,2,3 S 


sine [ | u £ x + v S y + w l z \ ] , 


(B.ll) 


At ? 

'4,5 = A T 


X, r = -nr- si ne[ | u ^ + v C y + w ^| + cU 2 + ^ + S z ) 1/2 ] • 


Consequently, there are five eigenvalues, and the time step is based on 
the maximum of them. 


A£ 

At- < ; 7 7 7 1/2 • 

1 |u c x + V 5 + w i z \ + CU 2 + 5y + C Z ) L/ 


(B.12a) 


Similar expressions can be found for At and At^ 

Ati 

At < 2 2 2 1/2 

11 |u ti + v n + w n 2 | + c(n x + n + n z ) 


(8.12b) 


At- < 7^ 7 7 7 T77 • 

C |u C x + V e + w C 2 | + c(C x + Cy + C 2 ) 


(B.12c) 
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B.2 Viscous Part 


where 


u = 


the 

viscous part, t 


5U 5F 


5t Fx 

P ^ 

/ 0 

pul 

1 

1 XX 

pv \ 

* F " \ " X xy 

pw 1 

/ T * 2 

pe) 

l V + ; 


0 

-T 


G = 


xy 


-T. 


yy 


-T 


yz 

q -d> 

M y V 


(B.13) 


0 

“ T > 

H = / -t 

“T 


XZ 


yz 

c 

zz 


q -<t> 
m 2 


Equation (B.13) can be transformed from physical coordinates to 
computational coordinates as 




H„ H 



(8.14) 


Similar to the previous case, Eq. (B.4) is split into three components. 
For example, the t component can be written as 


W +5 * F 5 + E y G 5 t5 z H !*°- 


After some algebraic manipulation, Eq. (B.15) can be written as 


(B.15) 


1$ + CA] Q rF + [B] q F _ + [C] q FF = 0 , 


FE 


XI 


Ch 


XI 


(B.16) 


where 
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[A] = 


0 

(4,+ti)^ y C x 

(4»+n)5 z 5 x 


._2 . 11 

U-n)S y +ng 

(4>+n)S z 5 y 


(4,+n)^ 

(«K»i)5 2 5 y 

U-lO^+WJ 11 

0 


Pg 


o 

o 

o 

o 

li 



0 

0 

0 

0 

0 


0 

{\-n)5 x Ti x +ng 12 

K xY v S\ 

<l>£> T] + |l£ TJ 
^x z ^z *x 

0 

II 

1 — 1 
CO 
LJ 

0 

sw, 

U-i*)? y i y +M 12 

^ y z 'y 

0 

j 

1 0 


•^V^Z 

U-n)S z T) x +ng 12 

0 


«g 


12 


Pg 


12 


j 

r 0 

0 

0 

0 

j 

0 

(X-H)^ x C x +Iig 13 


<\>Z c +us c 

XZ z X 

[C] H 

0 

«yWy 

(\-n)S y C y +ng 13 

<K y C 2 +n5 y C y 

1 

0 

«zWy 

4>S z +[il c 
z y z 

(\-^)c 2 c z +w 13 


aq 


13 


pg 


where 


X - 7 “ * 1*6 


* ■ “B ' 1 “ 


a 


y£ 

Pr p 2 (y-1) 


P ^ T r 1 

P Pr p(y-lj 
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Similar to the inviscid case, the following finite-Fourier series can be 


introduced as 


st '*1 
q = q • e • e 


ik,£ ik 9 n ik^C 


(B.17) 


Substitution of Eq. (B.17) into (B.16) results in the following 


[G] = [I - A - B - C] 


(B . 18) 


where 


2 6 1 

4At^ sin ig— 


A l 


T 


A, B = 


- sin 9j sin 0 2 At^ 




B , 


- sin 0. sin 0„ At_ 

c S Z~EZ ’ 


0J = k 1 M , © 2 = k 2 An , © 3 = k 3 AC . 

This system is stable if the largest eigenvalue is less than unity, this 
condition insures the decay of error. Therefore, the determinant of [G] 
can be set equal to zero as 



(B.19a) 


where 


G u 


= A. .+ 
TJ 


B u 


+ C, 


1J 
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MGgg-X) [(G 22 - ^) ^^ 33 “^) ( G 44 - X) ~ ^32~^ (^ 43 “^-)] 

- ^(632“^) (G^^-X) - (G 32 ~X) (G^2~X)] 

+ (G 24 -M [(G 32 -X) (G 43 -X) - (G 33 -\) (G 42 -X)]] = 0 . (B.19b) 

Equation (B.19b) can be solved for X as 


X 


1 


* 0 



(B.20) 


or 


x_ = p[- 


4 sin 


2 ”1 

T 11 


"ST 


sin 0j sin e 3 

_ 


12 , Sin 9 1 sin 0 3 13 

g + ^ g 


At. 

] W 


After all algebraic manipulations, the following are found 


At C < p P r (y- 1 ) 


yp [ 2 — - ’iLpl h + [ + ’V + ^2 "z, 


A? 


y y 

AC At) 


,c x c y + c v c v + c 7 Z 7I , 

+ I — AC AC — N • (B. 21 a) 


Similar expressions can be found for At and At„ 

T) C 


At < 


YP 


2 2 2 

T1 + T| +T) £ T] + £ *1 + £ Tl 

r 0 x 'y z ^ 1 x ’x ^y 'y ^ z 'z 

L c 5 + Tr TZ 


n p H r (Y-u L -^z - 1 AS An 


I^Y ^Y + % ^1/ + ^ , 1 

+ 1 irk — — U ■ < B - 21t >» 
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, t , YP r, C x ? + C y + «L ,<x\ + C y\* C z"z 

At C P P r (Y-l) 1 ^2 1 « 4 1 


<x 5 x + C y 'y + C z «z 

aT^C 


n 1 . <b 


Three more eigenvalues need to be computed from Eq. (B.21b) 
requires solution of a nonlinear equation for each point i 
computational space. 


.21c) 

which 
n the 
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